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In the study of gravitational waves (GWs), the stochastic background generated by 
compact binary systems are among the most important kinds of signals. The reason 
for such an importance has to do with their probable detection by the interferometric 
detectors [such as the Advanced LIGO (ALIGO) and Einstein Telescope (ET)] in the 
near future. In this paper we are concerned with, in particular, the stochastic background 
of GWs generated by double neutron star (DNS) systems in circular orbits during their 
periodic and quasi-periodic phases. Our aim here is to describe a new method to calculate 
such spectra, which is based on an analogy with a problem of Statistical Mechanics. 
Besides, an important characteristic of our method is to consider the time evolution of 
the orbital parameters. 

Keywords: Gravitational waves; stochastic background; double neutron stars. 

PACS Nos.: 04.30.Db, 95.85.Sz 


1. Introduction 

A big challenge of modern Astrophysics is the detection of gravitational waves 
(GWs), which would provide a new window to observe the universe. In the same 
way we use electromagnetic radiation such as X-rays, 7 -rays, infrared or visible light 
to study the astrophysical objects, we could well use the observation of gravitational 
radiation as an efficient tool in the search of the properties of the several kinds of 
astronomical objects. Specifically, we are concerned with the stochastic background 
of GWs generated by double neutron star (DNS) systems. The study of this kind of 
signal is important because it is among the most probable sources to be detected in 
the near future or, putting in other words, these signals could form a foreground for 
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the planned GW interferometers eLISA, BBO, DECIGO, Einstein Telescope (ET) 
and Advanced LIGO (ALIGO). 

Therefore, given the importance of the stochastic background generated by DNS, 
we use it as an example of application of a new method, described in the present 
paper, to calculate such spectra. The starting point is the equation^ 


^■BG — - - [ ^source ^ (1) 

Z'obs J 

where /ibg represents the dimensionless amplitude of the spectrum, i/ 0 b s is the 
observed frequency, dR is the differential rate of generation of gravitational radia tion 
and /iscmrce is the amplitude of the emitted radiation by a single source, namely^ ^ 



where /i is the reduced mass of the system, M is the total mass and g?l is the 
luminosity distance. 

It is worth mentioning that Eq. 0 was obtained from an energy-flux equation. 
This equation was first derived in a paper by de Araujo et al^ and was used in their 
various papers. In particular, in a subsequent papei^ they gave a more detailed 
derivation of this equation, showing its robustness. Also, although apparently simple 
it contains the correct and necessary ingredients to calculate the background of GWs 
for a given type of source. 

Given Eq. (JT]) and once h sou rce is known, we are focused on the calculation of 
the rate dR. First, we write this rate in the form 


dRdV , 
iR = dV dz 


( 3 ) 


where dV is the comoving volume element and z is the redshift. The element dV/dz 
is known from cosmology (see Subsection 12.11) and dR/d,V will be calculated by 
means of the new method considered in this paper. 

So, in order to describe the process to obtain the stochastic background, we 
organized this paper as follows: in Section [2] we show the population characteristics 
of DNS and the elements of cosmology necessary to the calculation of the spectra; 
in Section [3] we explain the method itself; the results are shown and discussed in 
Section [U and in Section [5] we present the conclusions and perspectives. 


2. Population Characteristics of DNS 

For the initial mass function (IMF), which is one of the ingredients of our calcula¬ 
tions, we adopt the Salpeter distribution, namely- 


4>(m) = Am 


( 4 ) 
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where A = 0.17 is the normalization constant and x = 1.7. The normalization of 
this IMF is obtained by means of 


rmf 

/ rruf>(m)dm = 1 (5) 

J mi 

where we are considering rrif = 125M© and rrii = O.IMq^I. Furthermore, we consider 
that neutron stars are generated by progenitors with masses ranging from 8Mq to 
25M©. It is worth pointing out that, although the choice of the values for these 
minimum and maximum masses of the progenitors are subject of disc ussio n, the 
range 8 — 25M© is usual in the literature (see, for example, Ferrari et alP-^l where 
the authors used these values). Further, discussions on this issue can also be found, 
for example, in SmarttP and Carrol & OstlicP. 

Besides the IMF, we need to adopt a star formation rate density (SFRD), which 
is another ingredient that appears in the calculation of the background. There are, 
in the literature, many alternatives for such ingredient. In particular, we adopt the 
SFRD derived by Springel & HernquislP^l, namely 


P*(z) 


0 a(z—Zm) 


— Pn 


/3 — a + ae@( z Zm '> 


( 6 ) 


where a = 3/5, /3 = 14/15, z m = 5.4 and with p m = 0.15M©yr 1 Mpc 3 fixing 
the normalization. It is worth mentioning that Eq. © was obtained considering a 
ACDM cosmology in a structure formation scenario, where the density parameters 
have the values Dm = 0.3, Db = 0.04 and Da = 0.7; where the subscripts M, B and 
A referes to matter, baryonic matter and cosmological constant, respectively. Be¬ 
sides, the authors used for Hubble’s constant the value of H 0 = lOOhkms _1 Mpc _1 
with h = 0.7 and it was considered a scale invariant power spectrum with index 
n = 1, normalized to the abundance of rich galaxy clusters at present day (as = 0.9). 

It is worth mentioning that one could argue that a different choice for the SFRD 
would modify significantly our results and conclusions. Since we are here mainly 
concerned with a new method to calculate the stochastic background of GWs, we 
leave the discussion concerning how different SFRDs affect the spectrum of the 
background of GWs, among other issues, to another paper to appear elsewhere. 

With the IMF and the SFRD at hand, it is possible to determine the formation 
rate of DNS. First, we consider the mass fraction A nsns that is converted into double 
neutron star systemic, namely 


Ansns — Pnsfp^ns (7) 

where (3 ns is the fraction of binary systems that survive to the second supernova 
event, f p gives us the fraction of massive binaries (that is, binary systems where both 
components could generate a supernova event) formed from the whole population 
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of stars and <f> ns is the mass fraction of neutron star progenitors that, in our case 
and using Eq. a, is given by 


p25 

^ns = / 4>(m)dm. (8) 

J 8 

Numerically, one obtains <h ns = 5.97 x 10 _3 Mq . Following a paper by Regimbau 
and de Freitas PachecdUl one has /3 ns = 0.024 and f p = 0.136. Using these results, 
the binary formation rate for DNS is given by 


s(^) — ^nsnsP* (z) • 


(9) 


Concerning the orbital parameters of the DNS, for the sake of simplicity, we 
consider a uniform period distribution given b}H^I 


f(p) = f (Pu - Pi)- 1 Pl<P<Pu 

1 1 0 otherwise 


( 10 ) 


where P u and P\ are the maximum and minimum periods, respectively. 

Since our main aim here is to present an alternative method to calculate stochas¬ 
tic background, the choice above is suitable for this purpose. 

However, in order to be used in the calculation of the spectra, Eq. should 
be rewritten in terms of the frequency. This is achieved by changing variables via 
f(P)dP = g(i'orb)eh' 0 rb, where the period P and the orbital frequency v or b are 
related to each other by P = Algebraic manipulations yield 


ff<>orb) = ( p 1 ) for P\<P<Pn (11) 

V 41 J V olb 

Now, note that the frequency z/ 0 rb undergoes time evolution. Following the paper 
by PeterJ33 one obtains 


^orb — 0 
Z7T 


(27rz/ orb ,o) 8/3 - -K(t - t 0 ) 


1 -3/8 


( 12 ) 


where K = 96miTO25c 5 G 5/,3 (mi+m2)~ 1 / 3 , mi and m 2 are the masses of the compo¬ 
nents of the system and t'orbfi is the initial frequency. Therefore, we should carry out 
a further change of variables in Eq. m in order to include the time dependence 
given by Eq. da. This is obtained by means of g(z' 0 rb,o)c^ / orb,o = JT(z' 0 rb)dz' 0 rb 
where z' or b,o is the initial frequency, which was associated with the variable j/ or b in 
CD and H(v or b) is the new distribution which, after some algebra, is written as 


orb) = ^ p 1 p t (^orb,o) 5/3 (l / orb) 11/3 


( 13 ) 
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Notice that it would be necessary to perform a further coordinate transformation 
in H (z'orb) in order to put it as a function of the emitted frequency v instead of 
t'orb- Such a transformation, given by v = 2t' or b, is trivial and all the equations will 
be written as functions of v from now on. Still concerning it is related to the 
observed frequency I'obs by means of v = u 0 \, a (l + z). 

2.1. Cosmology 

To perform the calculation of the background of GWs it is necessary to specify the 
cosmology and its corresponding parameters. Here we consider a flat universe. 

An essential quantity is the comoving volume element, which is given by 

^7 = 4tt (^) r 2 z F(n M ,^A,z), (14) 

where 


F(Qm, z) = — , =. 

-\/(l + z) 2 ( 1 + Hmz) — z( 2 + z)£Ia 

and the comoving distance r z for a flat universe reads 


r z 

r z = TT / F(fl M ,^A,z')dz' 
J o 


where the density parameters obey the following relations: 


(15) 


(16) 


12m = I2dm + and 12 m + 12a = 1 (17) 

The values of these parameters are those presented in the previous section. 


3. The method 

Since we are considering the time evolution of the orbital frequency of the systems, 
we need to take this issue into account in the derivation of the rate dR/dV in 
Eq. (I3|). We will see that the derivation of this rate comes down to count systems 
that reach a given frequency at a given moment of time. 

We derive this rate by means of an analogy with a problem of Statistical Me¬ 
chanics. In this problem, the aim is to calculate the number of particles that reach 
a given area A in a time interval dt , that is, the objective is to calculate the flux F 
of particles. Basically, this flux is calculated by summing all the particles inside the 
volume dV = Adx adjacent to the area A and that move towards A with velocity 
v = dx/dt , where v obeys a distribution function g(y) (see Fig.[T|). Hence, the sum 
is obtained by integrating over all the positive values of v. 

With some modifications, the method to calculate the flux F can be used to 
determine dR/dV. First, we substituted the spatial coordinate x by the frequency 
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Fig. 1. The volume element dV adjacent to A. 

v and the velocity v by the time variation of the frequency, which is defined by 
v v = dv/dt. So, the number of systems in the interval dv adjacent to a particular 
frequency v is given by 


V’O) 


J (p(v')dv’ 


(18) 


where (p(v) is the non-normalized distribution of frequencies. 

Considering that the distribution gives the number of systems which have 

v v in the interval dv v , the number of systems in dv and with values of v v in the 
interval dv v is given by 


^ = (jffiw) dv - (19) 

Now, the next step is to determine the forms of <^(f) and r](v u ). First, the distribu¬ 
tion <p[y) is written in the form 


<p{v) = J n nsns {t 0 )H(v)dt 0 (20) 

where to is the instant of birth of the systems, n nsns is the formation rate density 
of the DNS and H{y) is given by Eq. (fl3l) . 

In the deduction of Eq. m we consider initially from which we have 


dn = H(y)dv , 


( 21 ) 
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which is the fraction of systems originated at the time to an d that have frequencies 
in the interval dv. Now, using n nS ns) we can write explicitly 


dn 

dvdV dto 


— Bnsns{t 


Now, integrating over dto, we get 


( 22 ) 


dn 

dV 


^nsns (' to)H(y)dt 0 


dv, 


(23) 


where the expression in brackets is the number of systems per unit frequency in¬ 
terval and per comoving volume at given time (or redshift), which is the desired 
distribution function 

In Eq. (E3l) , the limits t and t min of the integral are related to the redshifts z and 
z SU p, respectively, with the usual expression found in any textbook on cosmology; 
further, ^ sup (see, e.g., Regimbau and Mandic 200fP31) is given by 


•^sup — 


if 


^obs ^ 


1 + ^max 


^max 

^obs 


otherwise. 


(24) 


It is worth mentioning that we are considering rwx = 900Hz (see next Section) and 

^max — 20. 

On the other hand, ri(v v ) will have a peculiar form. First, note that the derivation 
of Eq. (THU) yields 


CiV li 

v v = — oc v 3 (25) 

after some algebraic manipulations. Then, we conclude that there will be just one 
value of v v for each value of v , which allows us to write rj(vy) as a Dirac delta 
function: 

r)(v v ) = NS(vy - v„ tP ) (26) 

where N is the total number of systems and v U}P is the particular value of v„ 
corresponding to each frequency v. 

Now, noting that the denominator of the term between parenthesis in Eq. (usd 
is the total number of systems, using the function given by Eq. (12611 and changing 
the differential dv by means of the chain rule, Eq. m assumes the form 

dB = ( y I N6(v v - v VtP )dv v . (27) 


Integrating over v v and rearranging the result, we obtain 
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R = ( 28 ) 

where R is the number of systems per time interval dt. Recalling that the rate R is 
per comoving volume, we may write 


and Eq. © assumes the form 




dR 

dV 


(29) 


h 


2 

BG 


^obs 


h 


2 

source 


dR dV 
dVl + z’ 


(30) 


where we included the term (1 + z) in order to consider the time dilation due to the 
expansion of the universe. Further, using this amplitude we can obtain the spectral 
amplitude, which is given by: 

S h = ^. (31) 

Z'obs 


4. Results 

Recall that the main aim here is to see how the time evolution of the orbital fre¬ 
quency affect the spectra of the GW background. Then we considered two cases: 
first, we used 10 -6 and 10~ 4 Hz as the minimum and the maximum frequencies (cor¬ 
responding to periods of P u = 10 7 sec and Pi = 10 5 sec, respectively); in the second 
case we considered 10 -5 and 10“ 3 Hz (which give us P u = 10 6 sec and P\ = 10 4 sec). 
Note that here P u and P\ refer to the period of the waves, while in Eq. [13] such 
parameters refer to the orbital period. The transformation is trivial in this case. 

Figure m shows the spectra we generated using the distribution given by 
Eq. (1131) . Further, in order to point out the effects of the time evolution of the 
frequency on the results, we plotted in the same figure the curves corresponding to 
the case where the evolution is not considered. 

In both cases, we can observe the influence of the evolution of the frequencies 
on the form of the spectra. Particularly, if the evolution is taken into account, 
there will be a reduction in the amplitude of the spectra at the region of maximum 
initial frequency and a spread towards higher frequencies. Concerning the maximum 
frequ ency of the radiation emitted by the systems, we are assuming the value ^ max = 

900IlM 

Fig- 0 also shows for comparison the sensitivity curves for LISA, eLISA^^l, 
BBciS, decigcD!, ET^ and ALIGCPI The sensitivity curve for LISA may be 
found at http://www.srl.caltech.edu/^shane/sensitivity/, 

Although the spectra shown in Fig. © do not form foregrounds and therefore 
cannot be detected by single interferometric detectors, a suitable correlation of two 
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v(Hz) 

Fig. 2. Spectra generated with the use of Eq. {Hi. Note that in each case, in the regions of lower 
frequencies, there is superposition between the curves generated with (w.e.) and without evolution 
(n.e.). Also shown for comparison the sensitivity curves for LISA, eLISA, BBO, DECIGO, ET and 
ALIGO. 

or more of such detectors could, in principle, detect this background |20|21122 j n 
fact, the analysis of the detectability by cross-correlation, among other issues, are 
discussed in our other paper‘d to appear elsewhere. 

5. Conclusion and perspectives 

In this paper we shown an alternative method to calculate the background gener¬ 
ated by cosmological DNS during their periodic or quasi-periodic phases. We used 
an analogy with a problem of Statistical Mechanics in order to perform such a 
calculation, as well as taking into account the temporal variation of the orbital 
parameters of the systems. In this method we can easily change the distribution 
functions and the parameters without the need of modifying the formalism. 

We adopted here a plane period distribution to study the influence of the time 
evolution of the orbital frequency on the spectrum of GWs generated by DNS. 

In subsequent papers to appear elsewhere we will use the formalism developed 
here to: (a) calculate the background of GWs generated by black hole (BH) binaries 
and NS-BH binaries. Moreover, it will discussed, among other issues, how different 
SFRDs affect the spectrum of the background of GWs; (b) consider the GW spectra 
generated by compact systems in eccentric orbits; and (c) calculate the background 
of GWs generated by the coalescence of compact binary systems. In all these papers 
we will also discuss the detectability of the corresponding spectra by the present 
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and forthcoming GW detectors. 
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